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We generalise the Driven Diffusive Lattice Gas (DDLG) model by using a combination of 
Kawasaki and Glauber dynamics. We find via Monte Carlo simulations and perturbation stud- 
ies that the simplest possible generalisation of the equivalence of the canonical and grand-canonical 
ensembles, which holds in equilibrium, does not apply for this class of nonequilibrium systems. 
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^>^ . For statistical systems in thermodynamic equilibrium, the equivalence of different ensembles in the thermodynamic 

CT^ ' limit is a well-established result |l|] . Does this have an analogue in driven systems which display nonequilibrium steady 
0^ , states and transitions between them? Perhaps not, in general; but it is important to investigate when, if at all, such 
^~~^ ' an analogue might exist. One aspect of this problem has been studied in Ref. |g| in the context of a Gallavotti-Cohen- 
type symmetry in the large-deviation functional for driven stochastic systems such as a Driven Diffusive Lattice Gas 
(DDLG). We have studied another aspect of this problem in the context of a generalised Driven Diffusive Lattice 
Gas (DDLG), which is one of the simplest driven models in statistical mechanics with a transition between different 
O^ nonequilibrium steady states. We begin by recalling that the conventional DDLG (see below and Ref. Q) uses number- 
^^ , conserving (Kawasaki) dynamics [Q to update particle positions; i.e., it is the analogue of the canonical ensemble in 
I ■ equilibrium. To develop a grand-canonical analogue we generalise the DDLG to include a chemical potential fi and a 
^ , /I— dependent fraction of updates which use nonconserving (Glauber) dynamics H]; the remaining fraction of updates 
■^ use Kawasaki dynamics. We show the following: (1) our generalised DDLG is ideally suited to examining the simplest 

O/ , nonequilibrium analogue of the equivalence of canonical and grand-canonical ensembles; (2) even in this simple driven 
system, the canonical and grand-canonical ensembles are not equivalent. We arrive at this result by using Monte Carlo 
simulations to study our DDLG and perturbation theory to investigate a continuum version of it. We end with some 
qn^ ■ remarks about the relevance of our work to studies of phase coexistence in sheared mesogenic fluids |^ . 
0^ ■ It is useful to begin with a recapitulation of some elementary facts: The DDLG is based on a latttice-gas model in 

which the occupation variables Ui assume the values 1 or depending on whether a particle is present or not at the 
site i. Such a model is simply related to an Ising model [|6| defined in terms of the spin variables Si = {2ni — 1) by 
the Hamiltonian 
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where the exchange coupling J and the magnetic field H are related, respectively, to the pair potential V and the 
. !^ ■ chemical potential /i of the lattice gas, and {ij) are nearest-neighbour pairs of sites on a d— dimensional hypercubic 
K^ , lattice (we use a two-dimensional square lattice in our numerical studies). If J > 0, model (1) is ferromagnetic and 
Vh its lattice-gas analogue has an attractive interparticle interaction. The equilibrium phase diagram of model (1) is well 

known: In the temperature T and H plane there is a first-order phase boundary at H = along the line < T < Tc{d) 
which ends in a critical point at T = Tc{d); this first-order boundary shows up as a region of two- phase coexistence 
in a T — M phase diagram, where the magnetisation M is the Ising analogue of the lattice-gas density p; constant-M 
and constant-iJ ensembles are the analogues of the canonical and grand-canonical ensembles (we will use Ising-model 
and lattice-gas terminology interchangeably in this paper). As noted before, these ensembles are equivalent m and 
one can use standard thermodynamic relations to go from one to the other. In particular, to obtain the coexistence 
curve in the T — M phase diagram from the first-order boundary in the T — H phase diagram, we merely have to find 
the jump in the magnetisation across this phase boundary at all values of T < Tc{d). 

In the conventional DDLG, iJ = in model (1), the magnetisation M is fixed, since Kawasaki spin exchange 
is used in Monte Carlo updates, and a nonequilibrium steady state is maintained as follows: An "electric field" E 
is applied; this forces all particles (assumed identically charged) to move along its direction 1; periodic boundary 
conditions are used in this direction. In Monte Carlo simulations, one uses the algorithm of Metropolis et al M with a 
transition probability Min[l, e~''(^'^+''^)], where A7i is the change in energy because of the Kawasaki spin exchange 
and /3 = (fc^T)"^, with fc^ the Boltzmann constant. Note that the field E favours particles moving along its direction, 
disfavours the opposite, and does not affect jumps in transverse directions. Extensive studies [pi have shown that this 
DDLG exhibits two-phase coexistence for T < T^{d,E), where E = |E| and the superscript K stands for Kawasaki 
to help us to distinguish this critical temperature from the one we obtain below for our generalised DDLG. For the 



infinitely biased case, -B == oo, e.g., T^{d = 2,E = oo) ~ 1.35Tf (d = 2,E = 0), where T^{d = 2, £; = 0) is just 
the Onsager critical temperature for the two-dimensional Ising model in equilibrium. Critical exponents have been 
obtained for ii^ > ||] and one study ||] has investigated the coexistence curve in the T — M plane. 

We have generalised the DDLG by introducing Glauber spin-flip moves [|l^ in addition to the Kawasaki spin- 
exchange moves mentioned above. We choose the ratio Jq of the number of these Glauber moves to the total number 
of moves to be proportional to iJ^. Thus, as i/ ^ 0, /g ^ 0, in the simplest analytic way that is even in H . By virtue 
of these Glauber moves our generalised DDLG does not conserve the number of particles and thus provides a suitable 
extension of the grand-canonical ensemble for this nonequilibrium system. We might think naively that, as 7? — > 0, 
we regain the conventional DDLG with only Kawasaki updates. However, we must exercise caution here for there is 
some subtlety in the order in which limits are taken: since /g ~ H^ — > as if — > 0, we must run a Monte Carlo 
simulation for a time tss at least ~ H~^ so that the system experiences a large-enough number of Glauber moves 
and attains its true steady state; i.e., we must take the tss — > oo limit before we take iJ ^ (just as in equilibrium 
studies we take the thermodynamic limit (system size L — > oo) before we take the H ^ limit while calculating the 
magnetisation) . 

In our Monte Carlo simulations we use a square lattice of side L. In most of our studies E — oo and is applied along 
the +x direction. Thus jumps along this direction are always accepted, those in the —x direction are forbidden, and 
jumps along the ±y directions are not affected by E. We choose at random the spin that has to be updated, measure 
time in units of Monte Carlo steps per spin (MCS), and use random initial conditions. At each set of values of H and 
T we wait for the system to reach a statistical steady state, characterized by a steady mean value of the magnetisation 
per site {M{H, T) = {1/L'^) X]i ^i)j and then obtain data for average values of the quantities we measure. We obtain 
the coexistence curve in our dynamical grand-canonical ensemble by determining M{H,T) both as iJ t and H [Q, 
for T < T^^ , where the superscript GK indicates that this is the critical temperature for our generalised DDLG, 
which uses Glauber and Kawasaki spin updates. Curves of M versus T are shown for different values of H in Fig. 1. 
We use such curves to obtain the H 1 and H I limits of M{T, H) and thence the coexistence curve of Fig. lb (we 
show only the left half of this curve since it is symmetrical about M — Q or p — 1/2). 

Our coexistence curves for L = 16 and L = 32 (Fig. lb) are within error bars of each other, so finite-size corrections 
to our results are not significant except very near the critical point at p — 1/2, T = T^^ ~ 1.1. For comparison we 
have shown the coexistence curve obtained in Rcf. Q for the conventional DDLG, in which only Kawasaki updates 
are used; we also show the Onsager result for the two-dimensional Ising model in equilibrium. Figure lb illustrates 
two important features: (1) T^^ < T^ and the coexistence curve for our generalised DDLG is distinctly different 
from that for the conventional DDLG H; the former bows out to higher temperatures near T^^ , but then crosses the 
latter and subsequently lies below it. (2) The coexistence curve for our generalised DDLG is quite close to Onsager's 
result for the two-dimensional Ising model in equilibrium pj. We give a perturbative justification below. However, 
before we do this, it is useful to try to understand these results qualitatively. In our generalised DDLG, we approach 
the coexistence curve by taking the limits H 1 or H I 0. Thus, if T < T^^ , most spins assume the value sgn(iJ), 
and there are no macroscopically large interfaces as in the conventional DDLG. Consequently, the electric field E, 
which is the source of the nonequilibrium behaviour here, has a smaller effect in our generalised DDLG than it does 
in the conventional DDLG. This might well be the reason for the proximity of our coexistence curve to that of the 
two-dimensional Ising model in equilibrium. 
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FIG. 1. (a) The variation of the magnetisation M with the temperature T for our generalised DDLG with L — 32, E = oo, 

and the magnetic field H — ±0.05 (squares), H — ±0.02 (triangles), and H = ±0.01 (circles), (b) First-order phase-coexistence 
curves for the conventional DDLG (small diamonds) M, for our generalised DDLG (circles for L — 16 and triangles for L — 32), 
and for the two-dimensional Ising model in equilibrium (solid line). Note that p = (1 + M)/2. 



To obtain a more detailed understanding of our Monte-Carlo results, we have developed a Langevin or Time- 
Dependent Ginzburg-Landau (TDGL) model for our generalised DDLG. This is a simple extension of the Langevin 
model for the conventional DDLG |p[ ; since our purpose is merely to illustrate the phase-coexistence issues mentioned 
above, we restrict ourselves to a model in which all anisotropics, other than the driving electric field, are dropped. 
The Langevin equation for our model is 
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where A ~ (T — T^^) is negative in the ordered phase, and ■0 = + Mo, with M^ the mean-field magnetisation given 
by \Mo + (M/3!)Af^ — H . As in our lattice-gas model, the kinetic coefficients in this Langevin equation are such that 
the order parameter -0 is conserved if iJ = 0, and terms that lead to a violation of this conservation are proportional 
to H^ . The two noise terms 771 and 772 have zero mean and (77i(k, i)?7i(k', t')) = 2VoH^kBT5{k + 'k')5{t — t') and 
(?72(k, t)?72(k', t')) — 2k^kBT6(k + k')6{t — t'), where k and k' are wavevectors, k = |k| and t and t' are times. The 
variances of the noise terms are chosen such that, in the long-time limit, the Boltzmann distribution obtains ii E ~ 0. 
Since the current j^; produced by E must vanish if no holes or no particles are present locally, we make the simplest 
choice that satisfies these constraints, namely, j^; = (1 — ^^)E, which leads |^ to the term Edxi^^ in our Langevin 
equation with the spatial derivative along E (chosen to be parallel to the x axis). 

We now calculate M = (tp) perturbatively to 0{u,E'^) in the limits t -^ 00 and H ^ 0. The following diagrams 
contribute to this order: 






FIG. 2. Diagrams that contribute to the magnetisation M — {ip) to 0{u, E^) for model (2). Lines represent response 
functions and lines with a filled circle correlation functions. Vertices marked with a u have a factor of uH^ /2 associated with 
them; those with an E have a factor iE. 



Thus, to 0{u,E^) 



M = Mo + uMoksTai + uE'^Mo{kBTfa2 + uE^MoiksTfa^, 
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come from the loop integrals in Fig. 2; in order to compare with our lattice simulations we set the spatial dimension 
d = 2. 

The Langevin equation for the conserved case follows from Eq.(2) with H = 0. We set ip = (p + M, where, at the 
end of the calculation, we will find that M = Mo, to the lowest order in u. Hence 

^ ^ AV20 + ^V203 + ^m2v2,/) + ^Af V202 _ cWU + EMdx<t>'' + m- (5) 

The term 2EMdx4' has been eliminated by a Galilean shift. We calculate correlation functions by using the dynamic 
generating functional Q 

dtd'^x(f>V^(l) + i / dtd^x[[— - \W^(P -V2(/)^ - -M^V^cj) 
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- ^MW^<l>^ - cV''(/) -- Ed(t)^]$ + 0V20], (6) 

where (f> is the Martin-Siggia-Rose (MSR) conjugate variable ||l^. Order-parameter conservation implies that 4>{x,t) 
cannot respond to a spatially uniform magnetic field. In the two-phase regime phase separation proceeds via the 
formation of strips of up and down spins with the interfaces between these strips aligned, on average, parallel to E [gj. 
Thus the coupling to the field has the form / dtd'^x(l)V^h{x±) in the dynamical functional where subscript ± denotes 
the direction perpendicular to E. As in Ref. P] the equation of state is 
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where F is the vertex generating functional. A functional Taylor expansion of F about ipiq) = M5{q) yields 
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[0(92) + M5(g2)]...0(g„J<^(gi)...0feJxF"i"^(gi,...,(7„,,gi,...,g„J. (8) 

If we are only interested in the spontaneous magnetisation we set 

By retaining terms upto O(M^) we get 

_^[MF" + ^r- + ^F3i]=0 (10) 

which, when expanded to 0{u,E^), yields 



M = \/-^[l + ukeTbi + uE%2{kBTfbz + uE'^ikerfb^ + uE\kBTf[^i + A2 + A3)] 

Mo[l + ukgTbi + uEH2{kBTfh + uE^kBTfbs + uE\kBT)^{Ai + A2 + A3)], (11) 



where, as for the case of our generalised DDLG, the mean-field magnetisation Mo = y— 3!A74, and t*!, 62,^3 and A 
are the loop integrals in Fig. 3, To this order we can set M = Mo in the loop integrals. Hence we obtain 61 = oi, 
62 = 02, 63 = 03 and 



Ai 



A2 = 



1 



(PqdPqi 



1 



1 



qix{qx + qix) 



A J -2Aq2 + cg4 - 2X[q + qiY + c{q + qi)^ - 2Xqf 
1 /■ (fqd^qi 



4 



A3 = 



1 /■ d^^d^gi 



2Agf 



4 



-2A<72 + c<7^ -2A + cg2 -2A + cgj ' 

1 I gL 

-2Xq1 + cqf 2X + cq2 -2X + c{q + qi)^ ' 

1 1 qix{qx + qix) 



2XJ -2Aq2 + cg4-2A(9 + 9i)2 + c(g + gi)4-2Ag2 + cgf -2Ag2 + cqf -2A + cg2 -2A + eg? 



(12) 



Notice that sum of the diagrams contributing to F^^ and T'^^ to 0{uE^) vanish. This is a consequence of the invariance 
of our TDGL equations under r — > r — Ei with </) — > — 1/2. 

We now compare our TDGL results for the magnetisations of the generahsed DDLG and conserved cases. We find 
that there is an extra contribution from the last three diagrams Ai, A2,A3 in the latter; this is positive definite so 
I Mk I > I Mgk I • Of course if i? = both are the same as they must be by virtue of the equivalence of ensembles 
in equilibrium. Our analytical results agree qualitatively with our Monte Carlo results for 0.2 ^ p ^ 0.4 where 
the conventional DDLG coexistence curve lies above the one for our generalised DDLG (i.e., at a fixed value of 
T,pk > PGK or, equivalently, | Mk \>\ Mgk |); further away from this regime we must include higher-order terms 
in our functional Taylor expansion. In particular, we believe such terms are required to understand the crossing of 
the two coexistence curves in Fig. 1 for p < 0.2. Note also that quantitative agreement between our analytical and 
numerical results is not expected at criticality since our one-loop approximation can only yield mean-field exponents. 







FIG. 3. Diagrams contributing upto 0{u) and 0{uE^) to F^ 





FIG. 4. Diagrams contributing upto 0{u) to F^^ and F^^ 

In conclusion, then, we have shown that the simplest generalisations of grand-canonical and canonical ensembles 
are not equivalent for our generalised DDLG. Our study, though carried out on a very simple model, has important 
lessons for work on phase coexistence in systems such as sheared nematogenic fluids g] . Such studies have also found 
that constant-shear-rate and constant-stress ensembles yield different phase-coexistence boundaries. However, while 
determining such boundaries, the "chemical potentials" (defined as in equilibrium, i.e., as the derivative of a "free 
energy" with respect to particle density) in the two coexisting phases are equated. The lesson from our work is that 
this is valid only in the limit of very low shear rate (or E in our example); really we must equate dT /d<j) in the two 
coexisting phases; this will yield the equality of the chemical potentials in equilibrium but will have corrections at 
finite E (of 0{uE'^) to lowest order). 
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